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ABSTRACT 

We present three-dimensional Monte Carlo radiative transfer models of a very young (< 10^ 
years old) low mass (50 Mq) stellar cluster containing 23 stars and 27 brown dwarfs. The 
models use the density and the stellar mass distributions from the large-scale smoothed par- 
ticle hydrodynamics (SPH) simulation of the formation of a low-mass stellar cluster by Bate, 
Bonnell and Bromm. Using adaptive mesh refinement, the SPH density is mapped to the radia- 
tive transfer grid without loss of resolution. The temperature of the ISM and the circumstellar 
dust is computed using Lucy's Monte Carlo radiative equilibrium algorithm. Based on this 
temperature, we compute the spectral energy distributions of the whole cluster and the indi- 
vidual objects. We also compute simulated far-infrared Spitzer Space Telescope (SST) images 
(24, 70, and 160 /im bands) and construct colour-colour diagrams (near-infrared HKL and 
SST mid-infrared bands). The presence of accretion discs around the light sources influences 
the morphology of the dust temperature structure on a large scale (up to a several 10^ au). A 
considerable fraction of the interstellar dust is underheated compared to a model without the 
accretion discs because the radiation from the light sources is blocked/shadowed by the discs. 
The spectral energy distribution (SED) of the model cluster with accretion discs shows excess 
emission at A = 3-30 /im and A > 500 /im, compared to that without accretion discs. While 
the former is caused by the warm dust present in the discs, the latter is cause by the presence 
of the underheated (shadowed) dust. Our model with accretion discs around each object shows 
a similar distribution of spectral index (2.2-20 /im) values (i.e. Class O-III sources) as seen 
in the p Ophiuchus cloud. We confirm that the best diagnostics for identifying objects with 
accretion discs are mid-infrared (A = 3-10 /im) colours (e.g. SST IRAC bands) rather than 
HKL colours. 

Key words: radiative transfer - stars: formation - circumstellar matter - infrared: stars - 
brown dwarfs - accretion, accretion discs. 



1 INTRODUCTION 

Systematic investigations of young stellar objects (YSOs) in 
a star-forming cloud including comparative studies of theoreti- 
cal predictions and observations are important for understand- 
ing the stellar/substellar formation processes. Examples of well- 
studied star-forming clouds are the Orion Trapezium Cluster, 
NGC 2024, and the p Ophiuchus and Taurus- Auriga clouds. Re- 
cent observations of young clusters have been used to deter- 
mine circumstellar disc frequencies, disc mass distributions, ini- 
tial mass functions (IMFs) and evolutional stages of the ob- 
jects in the clusters. It has been shown that JHKL and HKL 
colour-colour diagrams are particularl y effective in identifying the 
presence of circumstellar discs (e.g. iKenvon & Hartmann 199?; 
McCaushrean et al. 1995; Ladaetal. 2000; Haisch et al. 2001). 
The spectral indices <Ladalll987ID or the slopes of observed spec- 
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tral energy distributions (SEDs) of YSOs from near- to mid- 
infrared wavelengths are used to classify SEDs, and the classifi- 
cation scheme is often related to the evolutionary s tages of YSOs 
(e.g. Adam s. Lada. & Shu 1987l lMvers et al. 1987). The distribu- 
tion of circumstellar disk masses in youn g clusters can be measured 
from millimetre continuum emission re.g. lEisner & CarDenteJ2003l 
for NGC 2024). Determining IMFs of young clusters requires evo- 
lution ary models (e.g. D'Antona & MazziteUi 1997; Baraffe et alJ 
1998) and infrared spectroscopy for spectral type identifications 
fe.g. Xuhman & Rieke^ 1999. for the p Ophiuchus cloud). 

IB ate. Bonnell. & Brom n] j2003h presented results from a very 
large three-dimensional (3-D) smoothed particle hydrodynamics 
(SPH) simulation of the collapse and fragmentation of a 50 Mq 
turbulent molecular cloud to form a stellar cluster. The calcula- 
tion resolved circumstellar discs down to 10 au in radius and 
binary stars as close as 1 au. Although some observational predic- 
tions, such as the IMF and binary fraction, may be gleaned directly 
from a hydrodynamical simulation of stellar cluster formation, the 
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principal observable characteristics (optical, near-IR, IR, and sub- 
millimetre images and spectra) require further detailed radiative 
transfer modelling. The density distribution of such hydrodynami- 
cal calculations is very complicated, and the corresponding radia- 
tive transfer must be also performed in full 3-D. 

There are two basic approaches to 3-D radiative trans- 
fer problems: grid based methods (e.g. finite differencing, 
short- and long- characteristic methods), and particle (pho- 
ton) based meth ods, i.e. Monte Carlo. E xampl es of the first 
kind are Stenholm . Storzer. & Wehrs a il99ll). [ Folini et alJ 
i2003h and ISteinacke^^^ mann. & He nnin d j20oj). Those' 
of th e second kind include Witt & Gordon n99u), "PasaniJ 
il99d),iWolf, Hennins, & Stecklum (1999), Harries (2000) and 



iKurosawa & H illiei (2001). The advantages of the second ap- 
proach are, for example, the flexibility to treat a complex density 
dis tribution and a coniplex scatte ring function. R e aders are referred 
to ISteinacker et al." i2003) and 'Pascucci et all l2004 ) for more 
extensive discussion on the advantages and disadvantages of these 
two different methods. 

The temperature of the interstellar and circumstellar dust in 
the cluster must be calculated in order to determine the source 
function of the dust emission. The radiative equilibrium temper- 
ature of the dust particles can be found using the Monte Carlo 
method (e .g. Lefevre, Bergeat, & Daniel 1982; Wolfetal. 1999; 
Lucy 1 19991: and iBiorkman & Woodll200lh . The technique used by 
Lucvl il999h takes into account the fractional photon absorptions 
between two events of a 'phot on packet' ; hence, it wo rks well even 
in the limit of low opacity. Biorkman & Wood ( 2001) used the im- 
mediate re-emission technique, in which radiative equilibrium is 
forced at each interaction with the dust. A photon is re-emitted 
immediately after an absorption event using a product of the dust 
opacity and the difference between the Planck function with a cur- 
rent temperature and that with a new temperature corresponding to 
the radiative equil ibrium of the dust that absorbed the photon. Un- 
like the method o f lLucvlJl999h . this does not require a temperature 
iteration if the opacity is independent of temperature. Alternatively, 
iNiccolini. Woitke. & Lopez (2003) used the ray-splitting method 
showing its effectiveness for both low and high optical depth me- 
dia. 

Here we aim to simultaneously resolve the dust on parsec 
and sub- stellar-radius spatial scales, whilst including multiple ra- 
diation sources. To overcome the resolution problem, we have 
implemented an adaptive mesh refinement (AMR) schem e in the 
grid production process of the TORU S radiative transfer ( Harries 
l200a). See alsolWolf et al]fl999) . Kurosawa & Hillier (2001) and 
eina ickeretalJJ200l^ for a similar griddi ng scheme in a radiative 
transfer problem. The method described by lLucvlJl999h is used to 
compute dust temperatures in our models. 

The objectives of this paper are: l.to compute observable 
quantities by solving the radiative transfer problem usi ng a com- 
plex d ensity distribution from the SPH calculation by iBate et all 
(l2003h : 2. to analyse the predicted obs ervational propert ies of the 
cluster generated in the simulation of Bate et al.' fTWJ) at a dis- 
tance of 140 pc which corresponds to the distance of nearby star- 
for ming regions such as Taurus-Auriga a nd the p Ophiuchus cloud 
(e.g.lBertout, Robichon, & Arenou|l999|). 

In Section 2, we describe the details of our models. The results 
of the model calculations are given in Section 3. The conclusions 
are summarised in Section 4. 
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Figure 1. Top: The scattering (dotted line), the absorption (dashed line), and 
the total (solid line) opacities of the grain model, described in Section IZSl 
are shown as functions of wavelength. Bottom: The corresponding albedo 
of the grains is shown as a function of wavelength. 



2 METHODS 

We calculate the main observable qua ntities of the low mass 
(50 Mq) cluster formation simulation by Bate et all i2003h . There 
are four basic steps in this process: 

(i) Create the AMR grid from the SPH particle positions and the 
density at each particle. The density is assigned to the grid cells 
during the grid construction. 

(ii) Calculate the luminosity and the effective temperature of 
each light source (in the SPH calculation) based on its mass using 
an evolutional model of you ng stars ( D ' Antona & Mazzitelli_ 1 998: 
ICensori & D'AntonJl99i) . 

(iii) Compute the temperature of the dust in the cluster using the 
Monte Carlo radiative equilibrium model of Lucv ( 1999). 

(iv) Compute the SEDs and the images using the Monte Carlo 
radiative transfer code, TORUS. 



2.1 The SPH model 



Bate et al. 1 2002 Eo02bL l2003h presented a numerical SPH sim- 
ulation of star cluster formation resolving the fragmentation pro- 
cess down to the opacity limit (^O.OOSM©). The initial condi- 
tions consist of a large-scale, turbulent molecular cloud with a 
mass of 50 Mq and a diameter of 0.375 pc (77400 au). The 
initial temperature of the cloud was 10 K, and its corresponding 
mean thermal Jeans mass was 1 Mq (i.e. the cloud contained 
50 thermal Jeans masses). The free-fall time of the cloud was 
tff = 6.0 X 10^^ s 1.90 X lOSr. Similar to the method 
used by "Ostrike r. Stone. & Gammid (|200?), they imposed an ini- 
tial supersonic turbulent velocity field on the cloud by generating a 
divergence-free random Gaussian velocity field with a power spec- 
trum P (k) = where k is the wave number. This was chosen 
to reproduce the observed Larson scaling relations ( Lars oi][T98ll) 
for molecular clouds. The total number of particles used in the sim- 
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ulations was 3.5 x 10^, making it one of the largest SPH calcula- 
tions ever performed. Approximately 95 000 CPU hours on the SGI 
Origin 3800 of the United Kingdom Astrophysical Fluids Facility 
(UKAFF) were spent on the calculation. 

We use the output from the final time step of the SPH calcula- 
tion as the input for the radiative transfer code. The time of the data 
dump is 0.27 Myr, at which the cluster contains 50 point masses 
(stars and brown dwarfs). Each SPH data point contains the posi- 
tion (x, y, z), the velocity components (Vx, Vy,Vz) and the density 
(p). The velocity information is not used in our calculation since 
we are only interested in continuum emission, but these data may 
be used in future calculations of molecular line emission. The total 
masses contained in the stars and the molecular gas are 6 Mq and 
44 Mq respectively. 



largest models split the grid on 27 levels, i.e. they incorporate a dy- 
namical range of 2^^ ^ 10^. The average optical depth across a cell 
is about 1.5 at 2 /xm which corresponds to the peak of a blackbody 
radiation curve (Bu) with the mean temperature (T = 2900 K) of 
the stars and brown dwarfs in our models. 



2.4 Temperature calculation 

The d ust temperature is computed using the method described by 
lLucvl (ll999). The material is assumed to be in local thermodynamic 
equilibrium (LTE) and in radiative equilibrium. The former indi- 
cates that the source function (Sx) is described solely by the Planck 
function, Bx, for all wavelengths A: 



Bx (T) 



(1) 



2.2 Source catalogue 

The SPH data provides the masses of the stellar objects (23 
stars and 27 brown dwarfs), ranging from 0.005 to 0.731 Mq 
(see Table [ij. Since the Monte Carlo radiative transfer code 
requires the luminosity (L), the effective temperature (Teff) and 
the radius (i?*) of each star, we computed them indirectly from 
evolutionary models. For a given age and mass of each star and 
brown dwarf, the luminosity and the tempe rature are interpolated 
from t he 1 998 updated versi on of data by D'Antona & Mazzitellil 
il998h and lCensori & D'Ant ona ( 1998) available on their website 
(http : / /www . mporzio .astro . it / "dantona/prems .html) 
Althou gh the actual age o f the stars and brown dwarfs in the SPH 
data of lBate et alJ J2003h ranges from -2 000 to -70000 years 
old, we make the pragmatic assumption that all the objects are 
0.25 Myr old because we beli eve that at young ag es the stellar 
radii are overestimated (c.f., iBaraffe et alJ l2002h . leading to 
unrealistically high luminosities. Adopting models with ages of 
0.25 Myr provides more plausible radii (and therefore luminosi- 
ties), while the change in temperature over this timespan is modest 
liaraffeet al. 2002). 

The radii are estimated from L = AnR^aT^^ where a is the 
Stefan-Boltzmann constant. The results are summarised in Table[l] 
The combined luminosity of all 50 sources is 16.8 Lq. 

We take the spectral energy distribution of the individual stars 
to be blackbodies. Although we recognise that the spectra of the ob- 
jects (particularly at the lowest masses) may be significantly struc- 
tured due to molecular opacity, the vast majority of the stellar flux 
is reprocessed by the circumstellar dust and the precise form of the 
input spectrum is unimportant. 



2.3 The grid construction (AMR) 

The algorithm use d to construct the gri d in th is paper is very sim- 
ilar to that used in lKurosawa & HillieJ i200lh . More detailed dis- 
cussion of the AMR grid construction and the data structure used 
in our model are given in Svmington (2004). Starting from a large 
cubic cell (with size d) which contains all the SPH particles, we 
first compute the density of the cell by averaging the density val- 
ues assigned to the particles. Then, the density is multiplied by the 
volume of the cell (d^) to find the total mass (M) in the cell. If the 
mass is larger than a threshold mass (Mth), which is a user defined 
parameter, this cell is split into 8 subcells with size d/2. If the mass 
of the cell is less than the threshold, it will not be subdivided. Re- 
cursively, the same procedure is applied to all the subcells until all 
the cells contain a mass less than the threshold (M < Mth). Our 



where T is the temperature of the dust. The latter indicates that the 
total energy absorbed by a volume element per unit time is exactly 
same as the total amount of energy emitted by the volume in unit 
time: 
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where kx and Jx are the thermal absorption coefficient and the 
specific mean intensity respectively. The expression can be de- 
rived by integrating the radiative transfer equation, dix/dl — 
pKx (Sx — Jx), over all wavelengths and solid angles (Q), us- 
ing the LTE c ondition (equation [I) and flux conservation (c.f., 
IChandrasekhaJ Tl 960 ) . By rewriting equation |2l using the Planck 
mean (hp) and the absorption mean (^ j), it becomes: 



where and Rj are defined as: 

^xB (T)^ d\ _^ KxBx (T) dX 
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and the (wavelength integrated) mean intensity J is 



J: 
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(3) 



(4) 



(5) 
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Using the A operator (e.g. 'Mihalas"l978'), the formal solution of 
the radiative transfer equation can be written as: 

J = A[S]=A[B (T)] . (7) 

This problem is often solved by iteratio n starting fro m some ini- 
tial temperature structure. As noted bv ILucvI i 19991) . equation |3l 
is indeed in the form of the formal solution; hence, we can use 
an iterative scheme to obtain the temperature. Note that using 
B (T) = aT^/iT in equation|3l we have: 



\crKp 



1/4 



(8) 



In addition, 'LucvI i 1 999|) explicitly used flux conservation as a con- 
straint when computing RjJ and in each iteration step. 
Our basic iteration scheme is summarised as follows: 

(i) Using the current temperature (Ti) from the i-th iteration 
step, estimate the va lues of Rj J and ^p for each cell using the 
Monte Carlo method iLucvll999l) . 
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Landscape table to go here. (Table 1 attached to the end of paper.) 

Table 1. 
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Figure 2. Top: the log of dust column density (in g cm~^) along x, y and z axes from left to right respectively (both Models I and II). Middle: The temperature 
on X = 0, y = 0, and z = planes from left to right (Model II). Bottom: The composite images of SST MIPS bands - 24 /im (blue) 70 /im (green) and 
160 /im (red) - along the x, y and z axes from left to right respectively (Model II). The black circles in the column density and temperature plots indicate the 
locations of the stars and brown dwarfs projected on each plane. 



(ii) Evaluate the new temperature (T^+i) using equation|8l 

(iii) Check for convergence. If the model has not converged, go 
back to (i). 



The model is considered to have converged when 



■AT, 



< 0.05 



(9) 



where AT, = Ti — Ti-i, and Ti is the mean temperature of all 
cells in the i-th iteration. 

Harries et al. (2004) presented test calculations for the dust 
temperature using our Monte Carlo model for a spherical shell 
(1-D) and a disc (2-D), and showed the results agree with 
those of well-estab lished existing codes Ilvezic et alj[l999l : 
iPascucci etal|200i). 



2.5 Dust model 

To calculate the dust scattering and absorption cross-sections as 
functions of w avelength, we have used the optical cons tants of 
IPraine & Led {1 984) for amorphous carbon grains and iHanneJ 
( 1988) for silicate grains. For simplicity, the model uses the stan- 
dard interstellar medium (ISM) power-law size distribution func- 
tion (e.g, Mathis. RumoL & Nordsieck_1977.) : 



n{a) (X a 



(10) 



for amin < a < amax with q — 3.5, amin = 5 x 10~ /xm and 
ttmin = 0.25 /xm. The relative number of each grain is assumed to 
be tha t of solar abundance, C/H~ 3.5 x 10~^ i Anders & Grevessd 
Il989l) and Si/H- 3.6 x 10"^ JGrevesse & Noelslll993l) which 
are similar to values found in the ISM model of iMathis et alJ 
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Figure 3. Log-scale temperature maps (in kelvins) from Model I (left) and Model II (right) on the 2; = plane. The circles indicate the positions of stars and 
brown dwarfs projected on the plane. The temperature of the dust shadowed by the accretion discs (left) is significantly lower than that of Model I. Although 
the same colour scale is used for both maps, the actual maximum dust temperatures for Model I and Model II are around 300 K and 1600 K respectively. 



J197I) and iKim. Mar tin. & Hendrvl ^199^). Similar abundances 
we re used in t he circumstellar disc models of ICotera et alJ J200lh 
and lWood et alj l2002b). Figure[l]shows the resulting opacities and 
albedo as functions of wavelength. 

Assuming spherical dust particles, the Mie- scattering phase 
matrix is pre-tabulated and used in our calculations, in- 
stead of using an ajjproxi mate scattering phase function (e.g. 
iHenvev & Greenstein 1941). We assume the same dust size dis- 
tribution in the ISM and the circumstellar discs, i.e. the circum- 
stellar discs do not contain larger grains. If larger grain sizes are 
chosen for the circumstellar dust, the wavelength dependency of 
the dust o pacity (Figu re Gl will be shallower than that of the 
ISM dust ( Wood et all 1 2002b). As a result, this may change the 
slope of model SEDs in sub-millimetre wave length range (c.f. 
iBeckwith et allToQclilBeckwith & SargenllQQll). 



2.7 Accretion disc model 



2.6 Images and SED calculation 

Once the temperature structure has converged, the source function 
is known everywhere in the model space (assuming LTE). Using 
this information, we simulate the observable quantities, namely the 
SEDs and the images. The Monte Carlo radiative transfer method is 
used again to propagate the photons, and project them on to an ob- 
server's plane on the sky. The basic method used here is presented 
in e.g.|Hillier ( 1991) and Harries (2000). 

The following sets of filters are used in our calculations: stan- 
dard /, H, K and L; SST Infrared Array Camera (IRAC) at 3.6, 4.5, 
5.8 and 8.0 /im; the SST Multiband Imaging Photometer (MIPS) at 
24, 70, and 160 fim. 



The SPH simulations o fiBateetaljliool do not resolve scale sizes 
less than ~ 10 au; hence, accretion discs very close to stars and 
brown dwarfs are not included. These missing accretion discs are 
potentially very important in the calculations of the SED and the 
images. Although large accretion discs exist around some objects 
in the SPH calculation, they do so only at distances > 10 au from 
the central object. To investigate the effect of the warm dust very 
close to the objects, we insert discs using the density described 
by the steady a-disc 'stand ard model' ( Shakura & Sunvaev 1973: 
iFrank, King, & Raine|2002|). 



Pd (r,z) = S(r) 



" ( 2h(r) ) 



2nh (r) 



(11) 



where r is the cylindrical radius expressed in units of the disc radius 
Rd, and h and z are the scale height and the distance from the disc 
plane, respectively. S is the surface density at the mid-plane. The 
mid-plane surface density and the scale height are given as: 

5Md 



S(r) 



-3/4 



where Md is the disc mass. 



h{r) = OmUdr 



9/8 



(12) 



(13) 



Note that the mid-plane surface density, predicted by the irradiated 
disk model of lD'Alessio et alJ (|^98t), has a slightly steeper radial 
dependency (S oc r~^). 

The inner radius of the disc is set to i?i = 6 which is 
the ass umed dust destruction radius of our model (^c.f. JWood et alJ 
l2002d) . The disc mass, Md, is assumed to be 1/100 of the cen- 
tral object's mass, and the disc radius (Rd) to be 10 au unless it 
has a binary companion. When an object is in a binary system. 
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Table 2. Basic Model Parameters 





Model I 


Model II 


Size of the largest cell 


2 X lO^s cm 


2 X IQis cm 


Size of the smallest cell 


1 au 


~ 0.1 Rq 


Discs within 10 au? 


No 


Yes 


Number of grid cells 


2,539,440 


25,223,192 


AMR subdivision levels 


16 


27 



the disc radius is assigned to be 1/3 of the binary separation (e.g. 
lArtv mowicz & Lubow" 1994). Finally, the orientations of the discs 
are assigned from the spin angular momentum of the objects in the 
SPH simulation which keeps track of the angular momentum of the 
gas accreted by the objects. 



3 RESULTS 

We present results from two different models: I. the density struc- 
ture is exactly the same as in the SPH simulation (i.e. there is no 
circumstellar material within lOau of any star) and II. the density 
structure is from the SPH simulation with the addition of accretion 
discs as described in Section ITTI The total number of cells used 
in Models I and II are 2,539,440 and 25,223,192 respectively. The 
corresponding cell subdivision levels are 16 and 27 for Models I 
and II, providing smallest cell sizes of 1 au and ^0.1 Rq respec- 
tively. Table|2l summarises the basic model parameters. 

3.1 Density and temperature maps 

We assume a gas-to-dust mass ratio of 100 fc.f. JLiseau et alll995h . 
and assign the dust density of each cell accordingly. The resulting 
dust column-density maps along x, y and z axes for Model II are 
given in the top row of Figure |2l On this scale, the column density 
maps for Model I look exactly the same except that the value of the 
maximum column density is lower. 

The results from the radiative equilibrium temperature calcu- 
lations for Models I and II are shown in Figure |3] Approximately 
500 CPU hours were spent on the SGI Origin 3800 of the UKAFF 
for computing temperatures for Model I, and 2400 CPU hours for 
Model II. Although the same colour scale is used for both models 
in Figure |3] the maximum dust temperatures found in the whole 
domain of the models are about 300 K and 1600 K for Models I 
and II, respectively. In both models, the average temperature of the 
ISM dust is around 20 K. The temperature maps for Model II on 
X =0, 2/ = 0, and z — planes are also given in the middle row 
of Figure |2| The angle averaged optical (at A = 5500 A) using 
600 random directions from the centre of the cluster to the outer 
boundary of the models is 93 {Av — 102). 

It is clear from the temperature maps of Model II that the 
presence of the circumstellar discs influences the morphology of 
the temperature structure. Discs affect the temperature of the dust 
even on a large scale (up to a several 10^ au). In Model II, 
the dust shadowed by the accretion disc has a lower temperature 
than that of Model I. The radiation field near the stars becomes 
anisotropic/bipolar when circumstellar discs are present. 

The density and the temperature maps around our typical star, 
object 4 (see Table [I), are shown in Figure |4] (upper right and left) 
for Model II. On the same scale, no structure will be seen in a sim- 
ilar plot for Model I as explained in Section 1231 The density map 



clearly shows that the cell size increases as the distance from the 
star (located at the centre) increases, and also as the distance from 
the mid-plane of the disc increases. In the temperature map, the 
hole around the star is created because the temperature of the low 
density material present in the background reaches 1600 K, the dust 
sublimination temperature; therefore, those cells are removed from 
the computational domain. We note that temperatures of the inner- 
most part of discs also can become greater than 1600 K, and the 
grid cells in this part of the discs will be removed. As a conse- 
quence, the inner radius may be slightly larger than an initial value 
{Ri = 6i?*). 

In the lower left of Figure |4] the density at the mid-plane of 
the disc is plotted as a function of the distance from the star. The 
slope of the density plot is ^ —15/8 as it should be according 
to equations [TT1IT31 Similarly, the temperature at the mid-plane is 
plotted as a function of the distance from the star in the lower right 
of the same figure. The slope at the outer part of the disc (0.3- 
10 au) is —0.45 which is consistent with that of Whitnev et all 
i2003bh . This is slightly smaller than the value for an optically 
thin disc (—0.4), with an absorption coefficient inv ersely propor- 
tional to wavelength, heated by stellar radiation rc.f.. lSDitzeJl978l: 
Kenvon et al. 1993; Chiang & Goldreich 1997). In the inner part 
of the disc (0.1-0.12 au), the slope quickly increases. A similarly 
rapid rise of the slope is seen i n the no-accretion-luminosity disc 
model of Whitnev et al J i2003d) as shown in their Fig. 6. 



3.2 Spectral energy distributions of the duster 

The SEDs from Models I and II are shown in Figure |5] along with 
the SED of the combined naked stars (the stars without the ISM and 
the circumstellar dust) and a 20K blackbody radiation spectrum. 
The observer is placed at a distance of 140 pc (on the +z axis) from 
the cluster. Both models show peaks around 2 /xm and 200 /xm. 
The first peak corresponds to that of stellar emission, and the latter 
corresponds to the emission from the relatively low temperature 
ISM dust. In Model I, the 200 /xm peak is prominent, indicating 
that the most of the ISM dust has T ^ 20 K. 

The most noticeable differences between the two models are 
the flux levels at A = 3-30 /im and A > 500 /xm. The excess emis- 
sion of Model II at A = 3-30 /xm is due to the warm dust present in 
the accretion discs. This emission arises from the reprocessing of 
photospheric emission by the inner discs into the observer's line-of- 
sight, and is dominated by a handful of objects (in particular object 
numbers 7 and 8, which constitute an ejected binary system). The 
excess emission in the far infrared (A > 500 /xm) in Model II is 
caused by the presence of a larger amount of cold (T < 20 K) dust. 
As we can see from Figure|3] a considerable fraction of the ISM in 
Model II is underheated compared to the dust in Model I because 
the radiation from the light sources is blocked/shadowed by the ac- 
cretion discs. The SED of Model II is very similar to the observed 
SED of a young star forming region (ISOSS J 20298-h3559-FIRl) 
with mass around 12OM0 presented by Krause et al. (2003). They 
showed that a Planck function with the temperature corresponding 
to the peak of the SED (T = 19 K) underestimates the flux levels 
observed at A < 60 /xm and A > 500 /xm. 

We note that the slope of the SEDs at sub-milimetere wave- 
lengths is sensitive to the wavelength dependency of the dust ab- 
sorption coefficient {n). If a shallower wavelength dependency of 
the dust opacity is introduced for the discs in Model II, the slope of 
the SEDs in the sub-millimetre wavelength range may change (e.g. 
iBeckwith et al|i;990;|Beckwith & Sargen|l99l|). 



8 R. Kurosawa et al. 







1000 



2 100 




Distance [AU] 



Distance [AU] 



Figure 4. The density (upper left) and the temperature (upper right) maps around a typical star (object 4 in Table [T], are shown for Model XL No material 
would be seen in a similar plot for Model I on the same scale due to the absence of small scale discs. In the temperature map, the hole around the star was 
created because the temperature of the low density material present in the background became more than 1600 K, the dust sublimination temperature, and the 
cells were removed from the computational domain. The density at the mid-plane of the disc is plotted as a function of the distance from the star (lower left). 
The temperature at the mid-plane is plotted as a function of the distance from the star (lower right). The temperature slope of the outer part of the disc (0.3-10 
au) is —0.45 (dash-dot line) which is consistent with that of Whitnev et al. (2003b). This is slightly smaller than the value for an optically thin disc (—0.4), 
with an absorption coefficient inversely proportional to wavelength, heated by stellar radiation ( Spitzei 1978; Kenvon et al. 1993; Chiang & Goldreiclj 
In the inner part of the disc (0.1-0.12 au), the slope quickly increases. A similar rapid raise of the slope is also seen in Fig. 6 of .Whitnev et aL i20 



3.3 Spectral energy distributions of stars and brown dwarfs 

Detailed st udies of is olated objects have b een performed by 
IWood etaljE)02a) andlWhitnev etlll J2003al). covering the evo- 
lution of the SED from the protostar through to the remnant disc 
stage. These models are an excellent point of comparison for our 
model SEDs, which are based on the density structures predicted 



for a 'realistic' star-forming cloud. The circumstellar material from 
the SPH calculation is however significantly more complex than the 
axisymmetric, idealized models cited above. 

To calculate the SEDs of individual objects, we have used 
the same density and temperature structures as in Section lTTI with 
some restrictions. Firstly, a cylinder of radius 50 au with its cen- 
tre passing through an observer (situated on the z-axis) and the 
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Figure 5. The spectral energy distributions (SEDs) from Model I (soHd 
line) & Model II (dash-dot line) compared with the 20 K blackbody ra- 
diation curve (dashed line) and the SED of naked (without the ISM and 
circumstellar dust) stars and brown dwarfs (dotted line). Both models show 
peaks around 2 iim and 200 fim. The first peak corresponds to that of stel- 
lar emission, and the latter corresponds to emission from the ISM dust with 
T ~ 20 K. The excess emission between 3 /xm and 30 /^m for Model II is 
caused by the warm dust in the accretion discs. The excess emission in the 
far infrared (A > 500 /im) seen in Model II is due to the emission from the 
increased amount of cold (T < 20 K) dust due to disc shadowing. 



centre of an object is considered. (The diameter of this cylinder 
roughly correspon ds to a 0.7 arcsec aperture at 140 pc). We note 
that W hitney etaP 12003 at) used 1000 au and 5000 au aperture sizes 
to compute colours of protostars, and found that the colours can 
change depending on the aperture size adopted. In our calculations, 
we use a smaller (50 au) aperture because the stellar densities in 
the cluster are such that larger apertures may include significant 
contributions from neighbouring objects. Secondly, the dust emis- 
sion, absorption and scattering outside of the cylinder are turned 
off. Thirdly, only the star under consideration can emit the pho- 
tons. With these restrictions, we have performed the Monte Carlo 
radiative transfer calculations for each object in Table [l] The re- 
sults are shown in Figure |6| Because of the restrictions used, some 
scattered flux contribution from the outside of the cylinder might 
be missing in the resulting SEDs, and the contamination due to the 
dust emission from the disc of a companion might be present in the 
SEDs if a star is in a binary system. Note that the binary pairs in 
the source catalogue (Table are 3-10, 7-8, 20-22, 44-42, 26- 
40, 39 -41 and 45-38. Readers are referred to table 2 of lBate etall 
i2003h for binary parameters (separations, eccentricities and so on). 
As we can see from the panels of Figure |6l the sources are deeply 
embedded in the cloud, and about half of them show little flux in 
th e optical. The SEDs are very similar to those of the Class model 
in'Whitnev et al.' ('2003 a). The silicate absorption feature at 10 /im 
becomes more prominent for objects with higher extinction. 

We have computed K band magnitude, the flux at 8 /am (Fg) 
and 24 /xm (F24) based on the SEDs, and extinction (Av)hy in- 
tegrating the opacities between each object and an observer. The 
results are placed in Table [l] For Model II, the inclinations of the 



discs with respect to an observer on +2; axis (that used in the cal- 
culation of the global SEDs) are also listed in the same table. The 
number of objects brighter than K = 20 are 26 and 30 objects for 
Model I and Model II respectively. 

Using the flux between 2.2-10.2 /xm, the SEDs in Fig- 
ure are classified based on the spectral index (a) as defined by 
Wilking, Lada, & Young ( 1989) : 

dlogjXFx) 
d log A 

We have used the monochromatic fluxes at K (2.2 /xm), L (3.6 /xm), 
M (4.8 /xm), N (10.2 /im) to compute the values of a via least- 
squares fits (e.g. Greene et al. 1994; Haisch et al. 2001). The clas- 
sification scheme of Greene et al. ( 1994) is adopted in our analy- 
sis. The sources with a > 0.3 are classified as 'Class I', —0.3 < 
a ^ 0.3 as 'flat spectrum', —1.6 < a ^ —0.3 as 'Class IF, and 
a < —1.6 as 'Class III' YS Os. In addition to this scheme, we clas- 
sify sources to be Class ( Andre et al.'l993') if the ratio of log XFx 
values at A = 160 /xm (SST MI PS) and at A = 2.2 /i m {K band) 
is greater than 3.0 (c.f.. Fig. 3 of lWhitnev et al 12003 ah . The results 
are placed in Table [l]for both models. Out of 50 objects, there are 
27 Class 0, no Class I, no flat spectrum, 1 1 Class II and 12 Class III 
objects for Model I, and there are 28 Class 0, 9 Class I, 4 flat spec- 
trum, 6 Class II, and 3 Class III objects for Model II. Note that even 
though all objects are < 0.07 Myrs old, there is a mixture of Class 
O-III objects. Thus, the class of an object does not necessarily re- 
late to its evolutionary stage. Whitnev et al. (2003a) also pointed 
out the degeneracy problem of the SEDs from Class I and Class II 
objects, and that their colours look very similar to each other for 
some combinations of disc inclinations and luminosities. 

Figure shows the distribution of the spectral index values. 
The upper row in the figure shows the histograms of the index val- 
ues from all 50 objects while the lower part shows those created 
from 'observable' (K < 20) objects. The histograms based on all 
objects peak around a = for Model I and II, and the distribu- 
tions are skewed. On the other hand, the number of the objects with 
K < 20 for Model II look more evenly distributed around a = 
(lower right in Figure 0. The latter is very similar to the spectral 
index distribution of the p Ophiuchus cloud presented in Fig. 3 of 
Greene et al. ( 1994), although we have to include fainter objects in 
order to have sufficient indices to make a comparison. 

3.4 Multi-band far-infrared images 

Figure |8] shows simulated images for the SST MIPS bands with 
central wavelengths (Ac) at 24, 70 and 160 /xm. These are idealised 
images in which the resolutions are limited only by the number of 
pixels used in simulations (500^ pixels), and are not degraded to 
the diffraction limit of the SST. The observer is placed at a dis- 
tance of 140 pc (on the -\-z axis) from the cluster. At this distance, 
the images subtend about 15.8 x 15.8 arcmin^. As predicted from 
the model SEDs (Figure|5}, most of the objects appear brighter for 
Model II in the 24 /xm image because of the warm dust emission 
from the circumstellar discs. Additional images for Model II with 
the observers on the +x and -\-y axis are also computed to create 3- 
colour, red (160 /xm), green (70 /xm) and blue (24 /xm), composite 
images. The results are shown in the third row of Figure |2l 

Although there is little effect of the circumstellar discs on the 
predicted SEDs around 70 and 160 /xm, the presence of the discs 
influences the morphology of the dust emission. In the lower part 
of the 70 and 160 /xm images of Model II, the 'butterfly' structure 
is caused by an almost edge-on disc. The structure resembles the 
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Figure 6. Atlas of individual spectral energy distributions (SEDs). The SEDs from Model I are shown in solid lines and those from Model II in dash-dot 
lines. The input source SEDs are presented in dotted lines. The numbers shown in upper right hand corner of each plot correspond to the object ID number in 
Table [l] To compute the SED from each star, the same density and temperature structures in section ITTI are used with some restrictions (see text for detail). 
As we can see from the SEDs, the sources are deeply embedded in a cloud, and about the half of them show little flux in the optical. The siHcate absorption 
feature at 10 fim becomes more prominent for objects with higher extinction. 



near infrared images of the ed ge-on discs around T Tauri st ars (e.g. 
HKTau C bv lStaDelfeldt~ l. 1998; HH 30 IRS by B urrows et alJ 
Il996h but while these butterflies are formed by scattering, those 
presented here are the result of thermal processes. Moreover, the 
scale size of the butterfly structure in the 70 and 160 /xm images 
of Model II is much larger (^ 10^ au) than that seen in the near 
infrared observations. The typical radius of the circumstellar discs 
around T Tauri stars ba sed on the near-infrared morphology is order 
of ^ 10^ au (e.g. see Burro ws et alJll994Euca s & Rochel ll998l : 
IStapelfeldt et aL .1998.) . Interestingly, a recent observation of the 



T Tauri star ASR 41 in NGC 1333 bv'Hodaop etall 12004 *) showed 
that the dark band in the reflection nebula around the star can be 
traced out to ^ 3000 au from the centre. Based on their radia- 
tive transfer models, they concluded that the large appearance of 
ASR 41 is probably caused by the shadow of a much smaller disc 
10^ au) being projected into the surrounding dusty cloud. 

To construct simulated images for actual observations by the 
SST, the images shown in Figure [S] are degraded to the diffraction 
limits of an 85 cm telescope; 7.6, 22.1 and 50.4 arcsec for 24, 70, 
and 160 /xm images respectively, by convolving them with a Gaus- 
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sian filter, and then the image pixels are binned-up to match the 
pixel scale of MIPS (2.55, 4.99 and 16.0 arcsec for 24, 70, and 
160 /im images respectively). The results are shown in Figure |9] 
The lower flux cut (the minimum value of the flux scale in each 
image) approximately corresponds to the MIPS limiting flux for 
each band (0.2, 0.5 and 0.1 MJy sr"^ for 24, 70 and 160 /im bands 
respectively). There is no significant change seen in the 24 /xm im- 
ages from the previous images while the 70, and 160 /im images 
are clearly degraded in resolution and sensitivity. No clear distinc- 
tion can be made between the two models (with or without the 
small scale discs) from the 70 /xm images. The major structures 
in the 160 /xm images appear to be the same in both models, but 
the emission from the isolated small structure in the lower half of 
the images is much weaker in Model II. According to the simulated 



24 /im images, 8 and 18 (out of 50) objects are detected above the 
flux limit for Models I and II, respectively. The simulated images 
at 70 and 160 /im show cloud structures with surface brightnesses 
that are 10-100 times the 557 detection limits. 

3.5 HKL and SST colour-colour diagrams 

Using the flux levels measured in the SEDs shown in Figure |6l 
we have constructed simulated HKL and SST IRAC (3.6, 4.5 and 
5.8 /im) colour-colour diagrams, and placed the results in Fig- 
ures [lO|and[n] We have also computed the intrinsic (de-reddened) 
colours of the objects in Model I and Model II by computing the 
SEDs without the ISM absorption (i.e. the opacity of the foreground 
dust is set to zero). The results are shown in the same figure for a 
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Figure 7. The distributions of the spectral indices for Model I (left) and Model II (right) computed from the flux levels at K, L, M and bands (2.2—10.2 jim) 
in the SEDs shown in Figure|6| While the upper plots include all 50 objects in our source catalogue (TableQ, the lower plots include only the objects with 
K < 20. The distribution for Model II with K < 20 (lower right) is very similar to that from the observation of the p Ophiuchus cloud bv Greene et alj 
J1994) shown in the inset. 



comparison. For the HKL colour-colour diagram, only the objects 
with K < 20 have been chosen, in order to simulate a deep pho- 
tometric observation. On the other hand, the objects with a flux 
greater than 10 /xJy (at 8 /xm) are used to construct the SST IRAC 
colour-colour diagram. 

As we can see from the left-hand plot in Figure fTol the ob- 
jects with circumstellar discs (Model II) are not well separated from 
those without discs (Model I), although the disc objects tend to be 
slightly redder than the discless objects. This can be understood 
from the SEDs presented in Figure |5] which shows no significant 
difference between the flux levels of Model I and Model II in the 
1-3 fim wavelength range. The intrinsic HKL colours are less scat- 



tered on the diagram. The disc obje cts seem to lie along the locus 
of cla ssical T Tauri stars given bv iMever. Calvet. & Hillenbrand 
<1997h : 

{H-K) = a{K-L) + b 

where a = 0.69 (±0.05) and b = -0.05 (±0.04). The least- 
squares fit of our disc object colours gives a = 0.51 (±0.27) and 
b = —0.11 (±0.37) which are in very good agreement with the lo- 
cus of Mever et al. ( 1997). This also indicates our disc SED models 
are reasonable. 

The disc objects and discless objects are well separated in the 
colour-colour diagram of the SST IRAC bands (in the left plot of 
Figure fm. Again, this can be roughly explained from the SEDs 
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Figure 8. Top: Idealised images for Model I at SST MIPS 24, 70, and 160 fim bands from left to right. Bottom: Same as the top images, but for Model 11. The 
units used for the colour scale in the figures are log of MJy sr~^. The cluster is placed at 140 pc from the observer, and the corresponding angular size of the 
images is 15.8 x 15.8 arcmin^. Most of the objects appear brighter in the Model II images due to the warm dust emission from the circumstellar discs. The 
presence of the disc influences the morphology of the dust emission in the 70 and 160 /im images. The butterfly structures caused by the almost edge-on disc 
can be seen in the lower part of the 70 and 160 fim images for Model II. 



in Figure |5l The SEDs from Model I and Model II start separat- 
ing from each other around A = 3 /im, and the difference between 
the flux levels increases as the wavelength increases until it reaches 
about 10 /xm. Moreover, the flux level of Model I decreases as a 
function of the wavelength and vice versa for that of Model II, in 
this wavelength range (3-10 /xm). The least- squares fit of the disc- 
less objects in the left plot of FigurelTTI gives: 

(loglO ^4.5 - logio ^3.6) = a (logio ^5.8 - log^o ^4.5) + b 

where a = 1.36 (±0.14) and b = 0.20 (±0.04). About 85 per cent 
of the disc objects are located to the right of this fit line. 

The intrinsic (de-reddened) colours of the objects are com- 
puted in the same way as was done for the HKL colours, and 
the results are placed in the right-hand plot of Figure [H] The 
least- squares fit of the disc objects is also shown in the same 
plot (and also in the left-hand plot). The slope and the inter- 
cept of the line are a = 0.78 (±0.50) and b = -0.04 (±0.05) 
respectively. While most of the disc objects are located above 
(log^Q F5.8 — logio ^4.5) = —0.29 and along the fit line, the disc- 
less objects are located below (logio ^s.s — logio ^4.5) = —0.29 
and along the same fit line. 

Observers (e.g . lAspin & Barsonvl Il994l: iLada et alJ l200Cl : 
iMuench et alJl200lh often use the JHK or HKL colour-colour di- 
agrams (either J-H vs H-K or H-K vs K-L) to find candidates for 



young stars or brown dwarfs with accretion discs. As we have seen 
in this section, a better diagnostic for identifying objects with ac- 
cretion discs is to use mid-infrared colours (e.g. SST IRAC bands) 
instead of the near-infrared colours. This may not be true for a more 
massive cluster since the massive stars could produce a larger frac- 
tion of high temperature dust. 



4 CONCLUSIONS 

We have presented three-dimensional Monte Carlo radiative trans- 
fer models of a very young low-mass stellar cluster with multiple 
light sources (23 stars and 27 brown dwarfs). The density struc- 
ture and the stellar distributions from the large-scale SPH simula- 
tion of Bate et al. (2003) were mapped onto our radiative transfer 
grid without loss of resolution using an AMR grid. The temper- 
ature of the ISM and the circumstellar dust wascomputed using 
the Monte Carlo radiative equilibrium method o f lLucvUl999h . The 
results have been used to compute the SEDs, the far-infrared SST 
(MIPS bands) images, and the colour-colour diagrams of this clus- 
ter. 

We find that the presence of circumstellar discs on scales less 
than 10 au (Model II) influences the morphology of the temperature 
structure of the cluster (Figure|3}, and can affect the temperature of 
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Figure 9. Top: Simulated ^S^rMIPS images for Model I at 24, 70, and 160 /im bands from left to right. Bottom: Same as the top images, but for Model 11. 
As in Figure[8| the colour scale values are in log of MJy sr~^. The images are degraded to the diffraction limits, 7.6, 22.1 and 50.4 arcsec for 24, 70, and 
160 ^im images respectively, by convolving them with a Gaussian filter. Further, the image pixels are combined to match the instrument resolutions of MIPS 
(2.55, 4.99 and 16.0 arcsec for 24, 70, and 160 iim images respectively). The lower flux cut (the minimum value of the flux scale in each image) approximately 
corresponds to the MIPS limiting fluxes (0.2, 0.5 and 0.1 MJy ~^ for 24, 70 and 160 /im bands respectively). No clear distinction can be made between the 
two models at 70 and 160 /xm. The dust emission from the circumstellar disc affects the number of observable stars in the 24 iim images. 
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Figure 10. The H-K vs. K-L colour-colour diagrams for the discless objects from Model I (stars) and the disc objects from Model II (circles). The left plot 
shows the colours of objects (with K < 20) computed from the SEDs given in Figure |6| and the right plot shows the intrinsic (de-reddened) colours of the 
same objects. The loci of intrinsic co lours for the main-sequence stars (solid line), the giants (dashed line), and the classical T Tauri stars (dash-dot line) from 
iMever. Calvet. & Hil lenbran? ('1997^ are also shown along with the reddening vectors (dotted lines) extending to the upper right from the edges of the loci. 
The least-squares fit of the intrinsic colours of the disc objects (Model II) are shown in thick solid line (right plot). 
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Figure 11. The simulated SST (IRAC) colour-colour diagrams of the discless objects from Model I (stars) and the disc objects from Model II (circles). The 
log^o of the fluxes at 3.6, 4.5 and 5.8 /im are used here instead of magnitudes. The left plot shows the colours of objects (with Fg > 10 fiJy) computed from 
the SEDs given in Figure|6| along with the least-squares fit line of discless objects (dashed line). About 85 per cent of the disc objects are located to the right 
of this line. The right plot shows the intrinsic (de-reddened) colours of the same objects along with the least-squares fit (solid line) of the disc objects. The fit 
(solid line) is also shown in the left plot. While most of the disc objects are located above logio ^5-8 — logio ^4.5 = —0.29 (dotted line) and along the fit 
line, the discless objects are located below logio ^5-8 — logio ^4.5 = —0.29 and along the same fit line. 



the dust on large scales (up to a several 10^ au). The dust shadowed 
by accretion discs has a lower temperature than the model without 
the discs (Model I). The radiation coming from within a few au of 
the light sources is anisotropic/bipolar because of the circumstellar 
discs. 

The cluster SEDs (Figure |5} from Models I and II both show 
peaks around 2 /im and 200 /xm. The first peak corresponds to that 
of stellar emission, and the latter corresponds to emission from ISM 
dust with T ^ 20 K. The excess emission of Model II between 
A = 3-30 /im is due to warm dust in the accretion discs. The excess 
emission in the far infrared (A > 500 /im) seen in Model II is 
caused by the colder (T < 20 K) dust. 

Assuming the distance to the cluster to be 140 pc, we have 
constructed simulated images (Figure |9j for a SST MIPS (24, 70, 
and 160 /im) observation using the appropriate diffraction limits, 
the sensitivity limits and the angular sizes of pixels. The emission 
at 24 /xm traces the locations of the stars and brown dwarfs very 
well. No clear distinction can be made between the 70 /im image 
from Model I and Model 11. The major structures in the 160 /xm 
images appear to be same in both Models, but the emission from 
the isolated small structure in the lower half of the images is much 
weaker in Model 11. In the simulated 24 /xm images, 8 out of 50 
objects are detected above the flux limit for Model I, and 18 out 
of 50 objects are detected for Model 11. We note that 15 of the 27 
brown dwarfs in the cluster (Model II) have i^T < 20 and would 
therefore detectable via deep imaging, while those that are too faint 
tend to be the youngest, most deeply embedded sources. 

Using the flux levels between 2.2 /im and 10 /im, the spec- 
tral indices of each SEDs were computed. The objects were then 
classified according to the spectral index values and the ratio of 
log AFa values at A = 160 /xm (SST MIPS) and at A = 2.2 /xm 
(K band). We have found that 54 per cent of the objects are classi- 
fied as Class 0, none as Class I, none as flat spectrum, 22 per cent 



as Class II and 24 per cent as Class III for Model L For Model 
II, 56 per cent of objects are classified as Class 0, 18 per cent as 
Class I, 8 per cent as flat spectrum, 12 per cent as Class II, and 6 
per cent as Class III. We also found that the spectral index distribu- 
tion of Model II ( with K < 20 objects) is very similar to that of 
the p Ophiuchus cloud observation by Greene et al. ( 1994j. Even 
though our objects are 0.07 Myr old, there are a mixture of Class 
O-III objects; hence, the class does not necessary relate to the ages 
of YSOs. 

According to the simulated HKL and mid-infrared (SST 
IRAC) colour-colour diagrams (Figures [iniandini, the disc ob- 
jects (Model II) tend to be slightly redder than discless objects 
(Model I) in the H-K vs K-L diagram, but the two populations are 
not well separated. This can be understood from the model clus- 
ter SEDs (Figure |5} that show no significant difference in the flux 
levels between A = 1-3 /xm. On the other hand, the disc objects 
and discless objects are clearly separated in the colour-colour di- 
agram of the SST IRAC bands (Figure [Til. As expected, we find 
that longer wavelengths are more efficient in detecting circumstel- 
lar discs. For example mid-infrared (A = 3-10 /xm) colours (e.g. 
SST IRAC bands) are superior to HKL colours. We find that the 
intrinsic colours of the disc objects (Figure [TH s how a distribution 
very similar to that found by Mever et al. ( 1997). 

The work presented in this paper can be extended to study how 
the observable quantities (e.g. colours of stars in a young cluster) 
evolves with time, by performing radiative transfer model calcula- 
tions with density structures from a hydrodynamics calculation at 
different times (ages). Growth of dust grain sizes in accretion discs 
may become i mportant in predic ting how the colours of objects 
evolve (e.g. seelHansen & Travisfl97llWood et all20023). 

One of the obvious next steps in modelling star formation is 
to combine SPH and radiative-transfer in order to more accurately 
predict the temperature of the dust. It is interesting to note that the 
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computational expense of the radiative equilibrium calculation per- 
formed here on one SPH time-slice (2400 CPU hours on UKAFF) 
is a significant fraction of that required to perform the complete 
SPH simulation (95 000 CPU hours on the same machine) which 
involves may thousands of time steps. Although it is not immedi- 
ately clear how often the temperature would have to be computed 
during a hydrodynamic simulation (the radiative-transfer calcula- 
tion would have its own Courant condition), it is apparent that a 
much more efficient method for calculating the radiation transfer is 
required before a combined SPH/radiative-transfer simulation be- 
comes tractable. 
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( * ) From the SPH calculation of Bate etal1l2003l . ^^^^^^^^^^ 

( * * ) Computed (at an age of 1/4 Myrs) from the isochrone data available on http : / /www .mporzio .as tro . it / "dantona/prems . html" See also lD'Antona & Mazzitelll^l99gt . 

(t) Estimated by using the L (column 3) and Tgff (column 4) in the Stefan-Boltzmann law. 

(I) Based on the SEDs of individual objects computed for an observer on -\-z axis at 140 pc. See section |33l 



